1 Import libraries

source("R/bioage_estimate_median.R")
source("R/gompertz_draw.R")
source("R/weibull_draw.R")
source("R/generate_data.R")
library(tidyverse)
library(MASS)
library(psbcGroup)
library(brms)

2 Generation of survival data

generate_data = function(pop_n, obs_n, p, g, ltname, seed, force_recalc, X_plots, a = exp(-9), b = 0.085, X_rho, β_rho, non_zero_betas, X_scale, active_hazard, betasep, bscale, betaplot = F) {
  
  print("Generating βs...")
  betas = generate_betas(p = p, g=g, rho = β_rho, rho_between = 0,
                         seed = seed, mu_u = betasep, mu_l = -betasep, beta_scale = bscale,
                         plot = betaplot, non_zero_groups = non_zero_betas, active_hazard = active_hazard)
  
  # betas$beta = betas$beta + 0.01
  
  
  print("Generating X...")
  X = generate_X(n = pop_n, p = p, g = g, rho = X_rho, seed = seed,
                 rho_between = 0,
                 X_plots = X_plots, scale = X_scale)
  print("Done!")

  
  
  print("Generating population lifetable...")
  lt = generate_population_lifetable(N_pop = pop_n, M = p, betas = betas$beta,
                                   X = X$X, filename = ltname,
                                   seed = seed, force_recalc = force_recalc, a=a, b=b)
  print("Done!")
  
  plot = ggplot(lt, aes(x = t, y = mrl)) +
    geom_line() +
    labs(
      title = "Population Lifetable MRL",
      x = "Chronological age (years)",
      y = "Mean Residual Life"
    ) +
    theme_minimal()
  print(plot)
  
  print("Generating data from lifetable...")
  df_sim = create_dataset(M = p, n_obs = obs_n, G = g,
                          gsize = (p/g), lt = lt, betas = betas$beta,
                          seednr = seed+1, X_plots = F, a=a, b=b, X_rho = X_rho,
                          X_scale = X_scale, followup = 20) # 20 year followup
  
  return(list(df = df_sim, true_betas = betas))
}

3 Functions to estimate MRL

3.1 Frequentist AFT

aft_reg = function(df_sim, p) {
  true_betas = df_sim$true_betas
  df_sim = df_sim$df
  

  
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]
  
  if (dim(df_sim_train)[1] < p) {
    return(NULL) # High dim - will not fit anyways!
  }
  
  # Create the formula for AFT regression
  predictor_vars = paste(c(glue("x{1:p}")), collapse = " + ")
  aft_formula = as.formula(paste("Surv(age_start, age_end, status) ~", predictor_vars))
  
  # Fit the AFT model
  fit_aft = tryCatch({
      fit_aft = eha::aftreg(
        formula = aft_formula, 
        data = df_sim_train, 
        dist = "gompertz",
        control = list(
          maxit = 100,
          trace = FALSE
        )
      )
      fit_aft
    }, 
    error = function(e) {
      message("AFT model error: ", e$message)
      NULL
  })
  
  # If model fitting failed, return NA
  if (is.null(fit_aft)) {
    return(NULL)
  }

  
  # get estimated parameters
  sigma_est <- exp(fit_aft$coefficients["log(scale)"]) # canonical parametrization
  tau_est <- exp(fit_aft$coefficients["log(shape)"]) # canonical parametrization
  a_est <- tau_est / sigma_est
  b_est <- 1 / sigma_est
  
  betas_est_aft <- fit_aft$coefficients[1:p]
  linpred_aft <- rowSums(sweep(df_sim_test[,1:p], 2, betas_est_aft, "*")) # SCALE
  
  for (i in 1:nrow(df_sim_test)){
    mrl_uncon <- integrate(gomp_baseline_surv, lower = (df_sim_test$age_start[i] * exp(linpred_aft[i])), 
                           upper=Inf, a = a_est, b = b_est)$value
    s_cond <-  gomp_baseline_surv(df_sim_test$age_start[i] * exp(linpred_aft[i]), a = a_est, b = b_est) 
    df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred_aft[i])
  }
  
  # Get the bias
  aftreg_coef_rmse = sqrt(mean((betas_est_aft - true_betas$beta)^2))
  
    
  # Get the RMSE
  aftreg_rmse = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  
  # MRL prediction plot
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for EHA AFT",
         subtitle = glue("RMSE = {aftreg_rmse}")) +
    theme_minimal()
    
  print(plt)
  
  # Plot difference between true and estimates betas
  true_betas$betahat = betas_est_aft
  
  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    # Vertical lines between pred and actual
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for AFT", subtitle = glue("RMSE = {aftreg_coef_rmse %>% round(4)}")) +
    theme_minimal()
  print(plt)
  
  return(list(prediction = aftreg_rmse, coefficients = aftreg_coef_rmse))
}

3.2 PSBC Bayesian AFT

psbc_reg = function(df_sim, p, g) {
  true_betas = df_sim$true_betas
  df_sim = df_sim$df
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]
  
  Y_train = cbind(df_sim_train$age_start, df_sim_train$age_end, df_sim_train$status)
  Y_test = cbind(df_sim_test$age_start, df_sim_test$age_end, df_sim_test$status)
  
  # Create covariate matrix X
  X_train = as.matrix(df_sim_train[, paste0("x", 1:p)])
  X_test = as.matrix(df_sim_test[, paste0("x", 1:p)])
  
  XC = NULL # CONFOUNDERS
  
  # Create group indicator for each variable
  grpInx = rep(1:g, each = p/g)
  
  # Set hyperparameters (PRIOR)?
  hyperParams = list(
    a.sigSq = 2, b.sigSq = 2,
    mu0 = 0, h0 = 10^6,
    v = 10^6
  )
  
  # Set starting values https://cran.r-project.org/web/packages/psbcGroup/psbcGroup.pdf
  startValues = list(
    w = log(Y_train[,1]),
    beta = runif(p, -0.5, 0.5), # uniform draws from realistic range
    tauSq = rep(0.1, p),
    mu = 0.1,
    sigSq = 0.1,
    lambdaSq = 0.1,
    betaC = rep(0.11, 0)
  )

  
  # MCMC parameters format __ CHANGE __ 
  mcmcParams = list(
    run = list(
      numReps = 5000,
      thin = 5, # Update params every n steps
      burninPerc = 0.2 # Discard the first ... steps
    ),
    tuning = list(
      mu.prop.var = 0.1,
      sigSq.prop.var = 0.005,
      L.beC = 50,
      M.beC = 1,
      eps.beC = 0.001,
      L.be = 100,
      M.be = 1,
      eps.be = 0.001
    )
  )
  
  fit_bayes_aft = aftGL_LT(Y_train, X_train, XC, grpInx, hyperParams, startValues, mcmcParams)
  beta_samples = fit_bayes_aft$beta.p

  # Convert to long format
  beta_long = beta_samples %>%
    as.data.frame() %>%
    mutate(iteration = row_number()) %>%
    pivot_longer(cols = -iteration, 
                 names_to = "beta_index", 
                 values_to = "value") %>%
    mutate(beta_index = factor(beta_index, levels = paste0("V", 1:ncol(beta_samples))))
  
  # Create boxplot
  plt = ggplot(beta_long, aes(x = beta_index, y = value)) +
    geom_boxplot(fill = "lightblue", alpha = 0.7) +
    geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
    labs(
      title = "Posterior Distributions of Beta Coefficients",
      x = "Beta Index",
      y = "Coefficient Value"
    ) +
    theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
  
  print(plt)

  
  # Extract betas using the maximum likelihood via density
  get_mode = function(x){
    xdist = density(x)
    mode = xdist$x[which.max(xdist$y)]
    return(mode)
  }
 # betas = fit_bayes_aft$beta.p %>% apply(MARGIN = 2, FUN = get_mode)
  betas = apply(fit_bayes_aft$beta.p, 2, mean)
  # Make linear predictor for MRL
  X_test = X_test %>% as.matrix
  linpred = X_test %*% betas
  
  # Change to lognormal
  # Create baseline lognormal survival
  # Implemented in R -> check log scales
  # dlnorm, use any of them 1 - F(x)
  # instead of gomp_baseline_surv no false, 1-
  
  
  # a = exp(-9)
  # b = 0.085
  
  # 1-plnorm(1:100, meanlog = fit_bayes_aft$mu.p[120], sdlog = fit_bayes_aft$sigSq.p[120])
  
  # Estimate mean from MCMC draws
  mu_est = mean(fit_bayes_aft$mu.p)
  sigSq_est = mean(fit_bayes_aft$sigSq.p)

  for (i in 1:nrow(df_sim_test)){
    mrl_uncon = integrate(function(x) 1 - plnorm(x, meanlog = mu_est, sdlog = sqrt(sigSq_est)), # Var to Sd via sqrt
                       lower = (df_sim_test$age_start[i] * exp(linpred[i])), 
                       upper = Inf)$value
    s_cond = 1 - plnorm(df_sim_test$age_start[i] * exp(linpred[i]), meanlog = mu_est, sdlog = sqrt(sigSq_est))
    df_sim_test$pred_mrl[i] <- (mrl_uncon / s_cond) * exp(-linpred[i])
  }
  
  # Get the RMSE
  pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  
  beta_RMSE =sqrt(mean((betas - true_betas$beta)^2))
  
  # MRL prediction plot
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for PSBC Bayesian AFT",
         subtitle = glue("RMSE = {pred_RMSE}")) +
    theme_minimal()
    
  print(plt)
  
  # Plot difference between true and estimates betas
  true_betas$betahat = betas
  
  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    # Vertical lines between pred and actual
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for Bayesian AFT", subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
    theme_minimal()
  print(plt)
  
  return(list(prediction = pred_RMSE, coefficients = beta_RMSE))

}

3.3 BRMS Bayesian AFT

brms_reg = function(df_sim, p, g) {
  # NO GROUPS
  true_betas = df_sim$true_betas
  df_sim = df_sim$df
  # df_sim$time_to_event = df_sim$age_end - df_sim$age_start
  
  # Train/test split 50/50
  train_indices = sample(1:nrow(df_sim), nrow(df_sim)/2)
  
  df_sim_train = df_sim[train_indices,]
  df_sim_test = df_sim[-train_indices,]

  # Create formula for brms
  preds = paste(paste0("x", 1:p), collapse = " + ")
  formula = as.formula(paste("age_end | cens(1 - status) + trunc(lb = age_start) ~", preds))

  # log(80)
  # > exp(log(80) - 1)
  # [1] 29.43036
  # > exp(log(80) + 1)
  # [1] 217.4625
  
  # Let's say that max should be 120...
  # exp(log(80) + 0.4) = 119.35, close enough

  
  priors = c(
    prior(horseshoe(df = 3, par_ratio = 0.1), class = b), # BETAS, horseshoe for shrinkage in high dim
    # https://rdrr.io/cran/brms/man/horseshoe.html, mainly refered to as alternative to spike and slab prior for shrinkage
    prior(normal(log(80), 0.4), class = Intercept), # about human life expectancy, 1 because i need a sd (uninformative?)
    prior(gamma(2,1), class = shape) #Hazard
  )


  # Fit Bayesian AFT with Weibull
  fit_bayes_aft = brm(
      formula = formula,
      data = df_sim_train,
      family = weibull(),
      prior = priors,
      warmup = 2000,iter = 4000, chains = 4, cores = 4
  )
  
  # Extract posterior for betas
  post = posterior_samples(fit_bayes_aft)
  
  # Extract beta coefs
  beta_cols = paste0("b_x", 1:p)
  beta_samples = post[, beta_cols]
  
  get_mode = function(x) {
    d = density(x)
    d$x[which.max(d$y)]
  }
  # betahat = apply(beta_samples, 2, get_mode)
  
  betahat = apply(beta_samples, 2, mean)
  betahat = -betahat # flip?
  
  # Set up linear predictor
  X_test = df_sim_test[,1:p] %>% as.matrix
  linpred = X_test %*% betahat
  
  # Calculate MRL using weibull
  shape_est = mean(post$shape)
  intercept_est = mean(post$b_Intercept)
  
  # Calculate MRL for each test observation
  df_sim_test$pred_mrl = NA
  
  for (i in 1:nrow(df_sim_test)) {
    # Weibull scale parameter for individual i
    scale_i = exp(intercept_est - linpred[i]) # "-" Due to flipping; this is personal hazard
    
    t_start = df_sim_test$age_start[i]
    
    # error handling for divergent integrals
    mrl_uncon = tryCatch({
      integrate(function(x) pweibull(x + t_start, 
                                   shape = shape_est, 
                                   scale = scale_i, 
                                   lower.tail = FALSE),
               lower = 0, upper = Inf)$value
    }, error = function(e) NA)
    
    if(is.na(mrl_uncon)) {
      df_sim_test$pred_mrl[i] = NA
    } else {
      s_cond = pweibull(t_start, shape = shape_est, scale = scale_i, lower.tail = FALSE)
      df_sim_test$pred_mrl[i] = mrl_uncon / s_cond
    }
  }
  
  
  # Calculate RMSEs
  pred_RMSE = sqrt(mean((df_sim_test$mrl - df_sim_test$pred_mrl)^2))
  beta_RMSE = sqrt(mean((betahat - true_betas$beta)^2))
  
  # MRL prediction plot
  plt = df_sim_test %>% ggplot(aes(x = mrl, y = pred_mrl)) +
    geom_point() +
    geom_abline(color = 'red') +
    labs(title = "True vs Predicted MRL for BRMS AFT",
         subtitle = glue("RMSE = {pred_RMSE}")) +
    theme_minimal()
    
  print(plt)
  
  # Beta comparison plot
  true_betas$betahat = betahat
  
  plt = true_betas %>%
    mutate(index = 1:p) %>%
    ggplot(aes(x = index)) +
    geom_segment(aes(xend = index, y = beta, yend = betahat), color = "gray") +
    geom_point(aes(y = beta, color = "True")) +
    geom_point(aes(y = betahat, color = "Predicted")) +
    scale_color_manual(values = c("True" = "black", "Predicted" = "skyblue"),
                       name = "Type") +
    ylab("β") +
    xlab("Index") +
    labs(title = "True vs Predicted β-coefficients for BRMS AFT", 
         subtitle = glue("RMSE = {beta_RMSE %>% round(4)}")) +
    theme_minimal()
  print(plt)
  
  return(list(prediction = pred_RMSE, coefficients = beta_RMSE))
}

3.4 Experiment function

run_experiment = function(pop_n = 20000, obs_n = 1000, p = 20, g = 4, ltname, seed, nzb, bsep, X_rho, plot_x = F, plot_betas = F) {
  
  cat("\n#### Data Generation\n\n")
  print("<<<<GENERATING SIMULATED DATA>>>>")
  df_sim = generate_data(pop_n = pop_n, obs_n = obs_n, p = p, g = g,
                       ltname = ltname, seed = 123,
                       force_recalc = F, X_plots = plot_x, a = exp(-9), b = 0.085,
                       X_rho = X_rho, β_rho = 0.9, non_zero_betas = nzb, X_scale = 0.005,
                       #active_hazard = 0.025, betasep = 2)
                       active_hazard = 0, # (0.05*sqrt(p))/nzb/bsep,
                       betasep = bsep, bscale = 1, betaplot = plot_betas)
  
  cat("\n#### Frequentist AFT Results\n\n")
  set.seed(seed+1)
  print("<<<<RUNNING AFT>>>>")
  rmse_aft = aft_reg(df_sim = df_sim, p = p)
  
  if (!is.null(rmse_aft)[1]) {
    cat("**AFT RMSE Results:**\n")
    cat("- Prediction RMSE:", round(rmse_aft$prediction, 4), "\n")
    cat("- Coefficient RMSE:", round(rmse_aft$coefficients, 4), "\n\n")
  } else {
    cat("**AFT Model failed to converge**\n\n")
  }
  
  cat("\n#### PSBC Bayesian AFT Results\n\n")
  set.seed(seed+1)
  print("<<<<RUNNING PSBC BAYESIAN>>>>")
  rmse_psbc = psbc_reg(df_sim = df_sim, p = p, g = g)
  
  if (!is.null(rmse_psbc)[1]) {
    cat("**PSBC RMSE Results:**\n")
    cat("- Prediction RMSE:", round(rmse_psbc$prediction, 4), "\n")
    cat("- Coefficient RMSE:", round(rmse_psbc$coefficients, 4), "\n\n")
  } else {
    cat("**PSBC Model failed to converge**\n\n")
  }
  
  cat("\n#### BRMS Bayesian AFT Results\n\n")
  set.seed(seed+1)
  print("<<<<RUNNING BRM>>>>")
  rmse_bms = brms_reg(df_sim = df_sim, p = p, g = g)
  
  if (!is.null(rmse_bms)[1]) {
    cat("**BRMS RMSE Results:**\n")
    cat("- Prediction RMSE:", round(rmse_bms$prediction, 4), "\n")
    cat("- Coefficient RMSE:", round(rmse_bms$coefficients, 4), "\n\n")
  } else {
    cat("**BRMS Model failed to converge**\n\n")
  }
  
  cat("\n#### Summary Comparison\n\n")
  
  results_df = data.frame(
    Method = c("AFT", "PSBC", "BRMS"),
    Prediction_RMSE = c(
      ifelse(is.null(rmse_aft), NA, rmse_aft$prediction),
      rmse_psbc$prediction,
      rmse_bms$prediction
    ),
    Coefficient_RMSE = c(
      ifelse(is.null(rmse_aft), NA, rmse_aft$coefficients),
      rmse_psbc$coefficients,
      rmse_bms$coefficients
    )
  )
  
  print(knitr::kable(results_df, digits = 4, 
                     caption = paste("Model Comparison for p =", p, ", n =", obs_n)))
  
  
  return(list(aft = rmse_aft, psbc = rmse_psbc, bms = rmse_bms))
}

4 Experiments

options(mc.cores = 4)

4.1 p = 20, g = 4

4.1.1 n = 250

rmses = run_experiment(p = 20, g = 4, ltname = "p20_may20_0", seed = 44, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 500, plot_x = T, plot_betas = T)

4.1.1.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: 1.32649924558738 Lower range of beta values pre-scaling: -1.65003671827946 Upper range of beta values pre-scaling: 3.73109001493681 Mean of betas post-scaling: -2.77555756156289e-17 Lower range of beta values post-scaling: -0.651799284773453 Upper range of beta values post-scaling: 0.374149793878827 [1] “Generating X…”

## Warning: `as.tibble()` was deprecated in tibble 2.0.0.
## ℹ Please use `as_tibble()` instead.
## ℹ The signature and semantics have changed, see `?as_tibble`.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: The `x` argument of `as_tibble.matrix()` must have unique column names if
## `.name_repair` is omitted as of tibble 2.0.0.
## ℹ Using compatibility `.name_repair`.
## ℹ The deprecated feature was likely used in the tibble package.
##   Please report the issue at <https://github.com/tidyverse/tibble/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

[1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 1985 Range of linpred values: -0.3918424 0.3903847

4.1.1.2 Frequentist AFT Results

[1] “<<<>>>” AFT RMSE Results: - Prediction RMSE: 3.78 - Coefficient RMSE: 0.1782

4.1.1.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Sun Jun 1 23:45:46 2025

iteration: 2000: Sun Jun 1 23:46:04 2025

iteration: 3000: Sun Jun 1 23:46:22 2025

iteration: 4000: Sun Jun 1 23:46:40 2025

iteration: 5000: Sun Jun 1 23:46:58 2025

PSBC RMSE Results: - Prediction RMSE: 13.3723 - Coefficient RMSE: 0.3162

4.1.1.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 56 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 10.1113 - Coefficient RMSE: 0.309

4.1.1.5 Summary Comparison

Model Comparison for p = 20 , n = 500
Method Prediction_RMSE Coefficient_RMSE
AFT 3.7800 0.1782
PSBC 13.3723 0.3162
BRMS 10.1113 0.3090

tion: 1 / 4000 [ 0%] (Warmup) Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 3.924 seconds (Warm-up) Chain 3: 2.168 seconds (Sampling) Chain 3: 6.092 seconds (Total) Chain 3: Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 3.993 seconds (Warm-up) Chain 4: 2.418 seconds (Sampling) Chain 4: 6.411 seconds (Total) Chain 4: Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 4.248 seconds (Warm-up) Chain 2: 2.386 seconds (Sampling) Chain 2: 6.634 seconds (Total) Chain 2: Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 4.322 seconds (Warm-up) Chain 1: 2.479 seconds (Sampling) Chain 1: 6.801 seconds (Total) Chain 1:

print(rmses)

4.1.2 n = 500

rmses = run_experiment(p = 20, g = 4, ltname = "p20_may20_1", seed = 44, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 1000, plot_betas = T)

4.1.2.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: 1.32649924558738 Lower range of beta values pre-scaling: -1.65003671827946 Upper range of beta values pre-scaling: 3.73109001493681 Mean of betas post-scaling: -2.77555756156289e-17 Lower range of beta values post-scaling: -0.651799284773453 Upper range of beta values post-scaling: 0.374149793878827 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 4022 Range of linpred values: -0.384812 0.3945516

4.1.2.2 Frequentist AFT Results

[1] “<<<>>>” AFT RMSE Results: - Prediction RMSE: 2.5561 - Coefficient RMSE: 0.1154

4.1.2.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Sun Jun 1 23:48:08 2025

iteration: 2000: Sun Jun 1 23:48:42 2025

iteration: 3000: Sun Jun 1 23:49:15 2025

iteration: 4000: Sun Jun 1 23:49:49 2025

iteration: 5000: Sun Jun 1 23:50:22 2025

PSBC RMSE Results: - Prediction RMSE: 14.0945 - Coefficient RMSE: 0.3623

4.1.2.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 9 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 4000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 2.7, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 9.0295 - Coefficient RMSE: 0.4226

4.1.2.5 Summary Comparison

Model Comparison for p = 20 , n = 1000
Method Prediction_RMSE Coefficient_RMSE
AFT 2.5561 0.1154
PSBC 14.0945 0.3623
BRMS 9.0295 0.4226

ion: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 13.207 seconds (Warm-up) Chain 2: 9.767 seconds (Sampling) Chain 2: 22.974 seconds (Total) Chain 2: Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 14.227 seconds (Warm-up) Chain 3: 10.23 seconds (Sampling) Chain 3: 24.457 seconds (Total) Chain 3: Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 399.354 seconds (Warm-up) Chain 1: 419.509 seconds (Sampling) Chain 1: 818.863 seconds (Total) Chain 1: Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 404.039 seconds (Warm-up) Chain 4: 419.185 seconds (Sampling) Chain 4: 823.224 seconds (Total) Chain 4:

print(rmses)

4.2 p = 60, g = 6

4.2.1 n = 250

rmses = run_experiment(p = 60, g = 6, ltname = "p60_may20_0", seed = 45, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 500, plot_x = T, plot_betas = T)

4.2.1.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: -0.741721151121596 Lower range of beta values pre-scaling: -5.83581464053189 Upper range of beta values pre-scaling: 3.66537305623651 Mean of betas post-scaling: 2.40548322002117e-17 Lower range of beta values post-scaling: -0.380425581789728 Upper range of beta values post-scaling: 0.384833071163438 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 2018 Range of linpred values: -0.3070134 0.3445294

4.2.1.2 Frequentist AFT Results

[1] “<<<>>>” AFT RMSE Results: - Prediction RMSE: 12.4042 - Coefficient RMSE: 0.3505

4.2.1.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Mon Jun 2 00:05:12 2025

iteration: 2000: Mon Jun 2 00:05:43 2025

iteration: 3000: Mon Jun 2 00:06:14 2025

iteration: 4000: Mon Jun 2 00:06:44 2025

iteration: 5000: Mon Jun 2 00:07:15 2025

PSBC RMSE Results: - Prediction RMSE: 13.1259 - Coefficient RMSE: 0.1826

4.2.1.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 675 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 2000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.68, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 15.6867 - Coefficient RMSE: 0.5214

4.2.1.5 Summary Comparison

Model Comparison for p = 60 , n = 500
Method Prediction_RMSE Coefficient_RMSE
AFT 12.4042 0.3505
PSBC 13.1259 0.1826
BRMS 15.6867 0.5214

ion: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 8.568 seconds (Warm-up) Chain 3: 3.318 seconds (Sampling) Chain 3: 11.886 seconds (Total) Chain 3: Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 9.881 seconds (Warm-up) Chain 2: 3.074 seconds (Sampling) Chain 2: 12.955 seconds (Total) Chain 2: Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 9.222 seconds (Warm-up) Chain 4: 6.785 seconds (Sampling) Chain 4: 16.007 seconds (Total) Chain 4: Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 214.456 seconds (Warm-up) Chain 1: 209.447 seconds (Sampling) Chain 1: 423.903 seconds (Total) Chain 1:

print(rmses)

4.2.2 n = 500

rmses = run_experiment(p = 60, g = 6, ltname = "p60_may20_1", seed = 45, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 1000, plot_betas = T)

4.2.2.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: -0.741721151121596 Lower range of beta values pre-scaling: -5.83581464053189 Upper range of beta values pre-scaling: 3.66537305623651 Mean of betas post-scaling: 2.40548322002117e-17 Lower range of beta values post-scaling: -0.380425581789728 Upper range of beta values post-scaling: 0.384833071163438 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 4039 Range of linpred values: -0.3365295 0.3515852

4.2.2.2 Frequentist AFT Results

[1] “<<<>>>” AFT RMSE Results: - Prediction RMSE: 7.0422 - Coefficient RMSE: 0.1871

4.2.2.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Mon Jun 2 00:15:55 2025

iteration: 2000: Mon Jun 2 00:16:52 2025

iteration: 3000: Mon Jun 2 00:17:50 2025

iteration: 4000: Mon Jun 2 00:18:47 2025

iteration: 5000: Mon Jun 2 00:19:44 2025

PSBC RMSE Results: - Prediction RMSE: 13.4759 - Coefficient RMSE: 0.1826

4.2.2.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 6000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 3.92, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 16.7679 - Coefficient RMSE: 0.5392

4.2.2.5 Summary Comparison

Model Comparison for p = 60 , n = 1000
Method Prediction_RMSE Coefficient_RMSE
AFT 7.0422 0.1871
PSBC 13.4759 0.1826
BRMS 16.7679 0.5392

on: 1 / 4000 [ 0%] (Warmup) Chain 4: Iteration: 1 / 4000 [ 0%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 13.193 seconds (Warm-up) Chain 4: 11.73 seconds (Sampling) Chain 4: 24.923 seconds (Total) Chain 4: Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 490.113 seconds (Warm-up) Chain 1: 526.625 seconds (Sampling) Chain 1: 1016.74 seconds (Total) Chain 1: Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 496.806 seconds (Warm-up) Chain 2: 525.894 seconds (Sampling) Chain 2: 1022.7 seconds (Total) Chain 2: Chain 3: Chain 3: Elapsed Time: 496.459 seconds (Warm-up) Chain 3: 526.366 seconds (Sampling) Chain 3: 1022.83 seconds (Total) Chain 3:

print(rmses)

4.3 p = 100, g = 10

4.3.1 n = 250

rmses = run_experiment(p = 100, g = 10, ltname = "p100_may20_0", seed = 46, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 500, plot_x = T, plot_betas = T)

4.3.1.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: 0.187187153779017 Lower range of beta values pre-scaling: -5.47642299999109 Upper range of beta values pre-scaling: 3.914142056722 Mean of betas post-scaling: -1.0547118733939e-17 Lower range of beta values post-scaling: -0.3662014611757 Upper range of beta values post-scaling: 0.236004069854284 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 1951 Range of linpred values: -0.2919669 0.3437821

4.3.1.2 Frequentist AFT Results

[1] “<<<>>>”

## AFT model error: Singular hessian; suspicious variable No. TRUE:
##  = 
## Try running with fixed shape

AFT Model failed to converge

4.3.1.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Mon Jun 2 00:38:25 2025

iteration: 2000: Mon Jun 2 00:39:14 2025

iteration: 3000: Mon Jun 2 00:40:02 2025

iteration: 4000: Mon Jun 2 00:40:51 2025

iteration: 5000: Mon Jun 2 00:41:39 2025

PSBC RMSE Results: - Prediction RMSE: 12.9332 - Coefficient RMSE: 0.1414

4.3.1.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 1 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 4000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 2.16, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 10.7303 - Coefficient RMSE: 0.1404

4.3.1.5 Summary Comparison

Model Comparison for p = 100 , n = 500
Method Prediction_RMSE Coefficient_RMSE
AFT NA NA
PSBC 12.9332 0.1414
BRMS 10.7303 0.1404

on: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 6.17 seconds (Warm-up) Chain 1: 5.43 seconds (Sampling) Chain 1: 11.6 seconds (Total) Chain 1: Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 285.684 seconds (Warm-up) Chain 2: 56.51 seconds (Sampling) Chain 2: 342.194 seconds (Total) Chain 2: Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 284.323 seconds (Warm-up) Chain 3: 254.485 seconds (Sampling) Chain 3: 538.808 seconds (Total) Chain 3: Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 290.607 seconds (Warm-up) Chain 4: 253.1 seconds (Sampling) Chain 4: 543.707 seconds (Total) Chain 4:

print(rmses)

4.3.2 n = 500

rmses = run_experiment(p = 100, g = 10, ltname = "p100_may20_1", seed = 46, nzb = 0.75, bsep = 4, X_rho = 0,
                       pop_n = 100000, obs_n = 1000, plot_betas = T)

4.3.2.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: 0.187187153779017 Lower range of beta values pre-scaling: -5.47642299999109 Upper range of beta values pre-scaling: 3.914142056722 Mean of betas post-scaling: -1.0547118733939e-17 Lower range of beta values post-scaling: -0.3662014611757 Upper range of beta values post-scaling: 0.236004069854284 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 4009 Range of linpred values: -0.3238358 0.3785866

4.3.2.2 Frequentist AFT Results

[1] “<<<>>>” AFT RMSE Results: - Prediction RMSE: 13.5896 - Coefficient RMSE: 0.2538

4.3.2.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Mon Jun 2 00:53:10 2025

iteration: 2000: Mon Jun 2 00:54:42 2025

iteration: 3000: Mon Jun 2 00:56:14 2025

iteration: 4000: Mon Jun 2 00:57:47 2025

iteration: 5000: Mon Jun 2 00:59:19 2025

PSBC RMSE Results: - Prediction RMSE: 13.694 - Coefficient RMSE: 0.1414

4.3.2.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 2 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 6.8806 - Coefficient RMSE: 0.1117

4.3.2.5 Summary Comparison

Model Comparison for p = 100 , n = 1000
Method Prediction_RMSE Coefficient_RMSE
AFT 13.5896 0.2538
PSBC 13.6940 0.1414
BRMS 6.8806 0.1117

on: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 17.762 seconds (Warm-up) Chain 3: 13.381 seconds (Sampling) Chain 3: 31.143 seconds (Total) Chain 3: Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 18.237 seconds (Warm-up) Chain 4: 13.371 seconds (Sampling) Chain 4: 31.608 seconds (Total) Chain 4: Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 18.496 seconds (Warm-up) Chain 1: 13.37 seconds (Sampling) Chain 1: 31.866 seconds (Total) Chain 1: Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 18.716 seconds (Warm-up) Chain 2: 13.34 seconds (Sampling) Chain 2: 32.056 seconds (Total) Chain 2:

print(rmses)

5 High dimensional example (p = n)

Since obs_n = 200, train test split of 50/50 means that it is trained on 100 observations.

rmses = run_experiment(p = 100, g = 4, ltname = "p60_may20", seed = 43, nzb = 0.75, bsep = 4, X_rho = 0, obs_n = 200)

5.0.0.1 Data Generation

[1] “<<<>>>” [1] “Generating βs…” Mean of betas pre-scaling: 1.27384106294797 Lower range of beta values pre-scaling: -2.12122062002442 Upper range of beta values pre-scaling: 4.07858990782366 Mean of betas post-scaling: 1.11022302462516e-18 Lower range of beta values post-scaling: -0.327775060966528 Upper range of beta values post-scaling: 0.204244829879737 [1] “Generating X…” [1] “Done!” [1] “Generating population lifetable…” [1] “Done!” [1] “Generating data from lifetable…” Nr of valid indices: 796 Range of linpred values: -0.3089132 0.3215815

5.0.0.2 Frequentist AFT Results

[1] “<<<>>>”

## AFT model error: Singular hessian; suspicious variable No. TRUE:
##  = 
## Try running with fixed shape

AFT Model failed to converge

5.0.0.3 PSBC Bayesian AFT Results

[1] “<<<>>>” iteration: 1000: Mon Jun 2 01:00:44 2025

iteration: 2000: Mon Jun 2 01:01:04 2025

iteration: 3000: Mon Jun 2 01:01:25 2025

iteration: 4000: Mon Jun 2 01:01:46 2025

iteration: 5000: Mon Jun 2 01:02:06 2025

PSBC RMSE Results: - Prediction RMSE: 12.6083 - Coefficient RMSE: 0.1414

5.0.0.4 BRMS Bayesian AFT Results

[1] “<<<>>>”

## Compiling Stan program...
## Trying to compile a simple C file

Running /Library/Frameworks/R.framework/Resources/bin/R CMD SHLIB foo.c using C compiler: ‘Apple clang version 17.0.0 (clang-1700.0.13.3)’ using SDK: ‘MacOSX15.4.sdk’ clang -arch arm64 -I”/Library/Frameworks/R.framework/Resources/include” -DNDEBUG -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/Rcpp/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/unsupported” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/BH/include” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/src/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppParallel/include/” -I”/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/rstan/include” -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include ‘/Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp’ -D_REENTRANT -DRCPP_PARALLEL_USE_TBB=1 -I/opt/R/arm64/include -fPIC -falign-functions=64 -Wall -g -O2 -c foo.c -o foo.o In file included from :1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Dense:1: In file included from /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/Core:19: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/library/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: ‘cmath’ file not found 679 | #include | ^~~~~~~ 1 error generated. make: *** [foo.o] Error 1

## Start sampling
## Warning: There were 2 divergent transitions after warmup. See
## https://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## to find out why this is a problem and how to eliminate them.
## Warning: There were 2000 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
## https://mc-stan.org/misc/warnings.html#maximum-treedepth-exceeded
## Warning: Examine the pairs() plot to diagnose sampling problems
## Warning: The largest R-hat is 1.57, indicating chains have not mixed.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#r-hat
## Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#bulk-ess
## Warning: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
## Running the chains for more iterations may help. See
## https://mc-stan.org/misc/warnings.html#tail-ess
## Warning: Method 'posterior_samples' is deprecated. Please see ?as_draws for
## recommended alternatives.

BRMS RMSE Results: - Prediction RMSE: 10.3288 - Coefficient RMSE: 0.1497

5.0.0.5 Summary Comparison

Model Comparison for p = 100 , n = 200
Method Prediction_RMSE Coefficient_RMSE
AFT NA NA
PSBC 12.6083 0.1414
BRMS 10.3288 0.1497

on: 1 / 4000 [ 0%] (Warmup) Chain 4: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 4: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 4: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 4: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 4: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 4: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 4: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 4: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 4: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 4: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 4: Iteration: 4000 / 4000 [100%] (Sampling) Chain 4: Chain 4: Elapsed Time: 2.536 seconds (Warm-up) Chain 4: 1.668 seconds (Sampling) Chain 4: 4.204 seconds (Total) Chain 4: Chain 2: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 3: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 1: Iteration: 400 / 4000 [ 10%] (Warmup) Chain 2: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 1: Iteration: 800 / 4000 [ 20%] (Warmup) Chain 3: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 3: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 2: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 1200 / 4000 [ 30%] (Warmup) Chain 2: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 2: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 3: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 3: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 3: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 3: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 3: Iteration: 4000 / 4000 [100%] (Sampling) Chain 3: Chain 3: Elapsed Time: 76.275 seconds (Warm-up) Chain 3: 1.773 seconds (Sampling) Chain 3: 78.048 seconds (Total) Chain 3: Chain 2: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 2: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 2: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 2: Iteration: 4000 / 4000 [100%] (Sampling) Chain 2: Chain 2: Elapsed Time: 74.53 seconds (Warm-up) Chain 2: 12.725 seconds (Sampling) Chain 2: 87.255 seconds (Total) Chain 2: Chain 1: Iteration: 1600 / 4000 [ 40%] (Warmup) Chain 1: Iteration: 2000 / 4000 [ 50%] (Warmup) Chain 1: Iteration: 2001 / 4000 [ 50%] (Sampling) Chain 1: Iteration: 2400 / 4000 [ 60%] (Sampling) Chain 1: Iteration: 2800 / 4000 [ 70%] (Sampling) Chain 1: Iteration: 3200 / 4000 [ 80%] (Sampling) Chain 1: Iteration: 3600 / 4000 [ 90%] (Sampling) Chain 1: Iteration: 4000 / 4000 [100%] (Sampling) Chain 1: Chain 1: Elapsed Time: 118.518 seconds (Warm-up) Chain 1: 104.365 seconds (Sampling) Chain 1: 222.883 seconds (Total) Chain 1:

print(rmses)

Small tests:

6 Monte carlo (TODO)

# Parameters grid:
# p
# g
# 
# MC = replicate(10000, {
#   
# })